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Abstract 

A numerical approach to Ginzburg-Landau (GL) theory is demonstrated and we review its apphcations to several 
examples of current interest in the research on superconductivity. This analysis also shows the applicability of 
the two-dimensional approach to thin superconductors and the re-defined effective GL parameter k. For two-gap 
superconductors, the conveniently written GL equations directly show that the magnetic behavior of the sample 
depends not just on the GL parameter of two bands, but also on the ratio of respective coherence lengths. 
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Whenever a new scientific discovery is made, re- 
searchers must strive to explain it theoretically. In 
the case of superconductivity, it took more than two 
decades after its experimental discovery before the 
London theory was developed [1]. However, the Lon- 
don theory treats vortices as point-like objects and 
does not take into account the finite size and the 
inner structure of the vortex. In 1950, Landau and 
Ginzburg [2] developed a phenomenological theory, 
which combined Landau's theory of second-order 
phase transitions with a Schrodinger-like wave equa- 
tion. Over the past 50 years, this theory had great 
success in explaining macroscopic properties of su- 
perconductors (such as the division of superconduc- 
tors into two categories now referred to as type-I 
and type-II, and the description of the mixed state 
of type-II superconductors [3]), but it was also im- 
mensely useful for description of mesoscopic super- 
conducting samples [4,5]. Even the full, microscopic 
Bardeen, Cooper, and Schrieffer (BCS) theory of 
superconductivity reduces to the Ginzburg-Landau 
(GL) theory close to the critical temperature [6]. 



The GL theory extends Landau's theory of 
second-order phase transitions [7] to a spatially 
varying complex order parameter ■0(^) which is 
nonzero at T < Tc and vanishes at T > Tc through 
a second order phase transition. The resulting gra- 
dient term is made gauge-invariant by combining it 
with the vector potential A{r) where V x A{r) = 
h{r) is the local magnetic field. 

The two Ginzburg-Landau equations are obtained 
by minimization of the GL free energy functional 
J^IV", A} with respect to ?/> and A 



K'^{h-Ha)^ dV, 



(1) 
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where n is the GL parameter given as a ratio of mag- 
netic penetration depth A and the coherence length 
and Hq denotes the applied magnetic field. Eq. 
(1) is given in dimensionless form, where all dis- 
tances are measured in ^(T), the vector potential A 
in cft/2e^, the magnetic field H in Hc2 = ch/2e£,^, 
and the order parameter ip in i/jq = y'— a//?, with 
a, (3 being the material dependent coefficients. 
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Every part of Eq. (1) describes some physical 
property. In principle, it is possible to introduce 
some extra terms in the energy functional in order 
to describe the superconducting state deeper in the 
superconducting phase (sec, for example, Ref. [8]), 
but the achieved corrections arc very small and are 
rarely considered. The first part of Eq. (1) is the 
expansion of the energy difference between the su- 
perconducting and normal state for a homogeneous 
superconductor in the absence of an applied mag- 
netic field ^ near the zero-field critical temperature 
Tco- The coefficient a is negative and changes sign as 
temperature is increased over Tcq {a cc (T ~ Tco))) 
while [3 is a positive constant, independent of tem- 
perature. Therefore, the Cooper-pair density corre- 
sponding to temperatures below Tco and in absence 
of magnetic field is IV'oP — ~ct/f3. 

The next term in Eq. (f ) is clearly the kinetic en- 
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1 

2m* 



-ifiV ■ 



ergy of the Cooper-pairs 
where m* is the effective mass of a Cooper-pair. It 
describes the energy cost when the superconducting 
density is non-homogeneous. 

The last term in Eq. (1) describes the energy of the 
magnetic field of the supercurrents, which measures 
the response of the superconductor to an external 
field and is nothing else than the difference between 
the local and applied magnetic field. Note that for a 
superconductor in an external field, the equilibrium 
state is not defined by the Helmholtz free energy but 
the Gibbs free energy. 

The phenomenological GL theory is one of the 
most elegant and powerful concepts in physics, 
which was applied not only to superconductivity 
(see textbooks [9-12]) but also to other phase transi- 
tions, to nonlinear dynamics, to dissipative systems 
with self-organizing pattern formation, and even to 
cosmology. In what follows, we describe a numer- 
ical approach towards solving the GL equations, 
with several particular modifications for special 
problems in superconductivity. 



Numerical approach in general 

In what follows, we describe a numerical method 
used to solve the Ginzburg-Landau equations, whose 
solution minimizes the free energy. Using the dimen- 
sionless variables as explained above and the Lon- 



^ It can be shown from the microscopic theory that only 
even powers of appear in this expansion. 



don gauge, divA = 0, GL equations can be written 
in the following form: 

(-iv-A)%=v(i-i^r), (2) 

-K^AA= -^{^*\/i;-^py^*)-\'il;f A. (3) 

The Neumann boundary condition on the sample 
surfaces takes the form 



boundary 



= 0. 



(4) 



For the fixed applied magnetic field, we solve 
the two coupled Ginzburg-Landau equations self- 
consistently. Both equations are solved using the 
link variable approach [13] for a finite-difference 
representation of the order parameter and the vec- 
tor potential on a uniform cartesian space grid 
{x, y) with a typical grid spacing of less than 0.1^. 

For a given applied magnetic field, we start from 
the applied vector potential as initial condition in 
our calculation, as if no superconductor is present. 
The first step is to solve the first GL equation (2). 
According to Kato et al. [13], Eq. (2) can be written 
as 

_ 

dt ~ ' 



A] V + 



(l^l'-l)^ 



+ /(r,i),(5) 



where time relaxation is included on the left side, 
and /(r, t) is a dimensionless random force. It is es- 
sential to put the gauge field A on the links of the 
computational lattice, which is achieved by intro- 
ducing the link variables between r\ and r2 as 



C/;^''-^ = exp 



r2 

i j A/n (r) • dp. 

L fi 



(6) 



with /J = x,y, z. 

In our calculation, the whole system is mapped 

on a rectangular grid. The first term from Eq. (5) is 
discretized as (the index j denotes the lattice point 
of interest) 



1 



(7) 



+Ali^j+iA^V^ijj = — (-2iA^f/^V^Vi 

-iU^,M^,A,-iAl) + U^S/l^,). 

After substituting V^C/^ = -iAfj,W^ and V^?/^ = 
—iU^{\/ fiAfj^ — iA'jJ, and some trivial transforma- 
tions, we obtain 



2 



(8) 



Finally, for /i = a; (analogously for ij. = y,z) we have 



i 



-A 



(9) 



The discretized Ginzburg-Landau equation can be 
now written in full as 



at 



+ 



(10) 



where different indices denote adjacent grid points 
to point j along three axis. In general, this approach 
works for ax ^ ay i= a^- 

Using Eq. (10), with chosen initial vector poten- 
tial, we solve for the value of the order parameter ^ 
in every grid point. These values we can implement 
to calculate the local current densities jx,y,z- The 
right side of Eq. (3) can be written as 



^=2 



,(11) 



where again link variable approach comes into play 
through similar transformations as in Eqs. (7)-(10) 



:Vx - Ax tpj 



-»77jV.(C/iV,) 
Ui 



(12) 



From the supercurrents a new value for the vec- 
tor potential can be calculated using the second GL 
equation. This is a Poisson-type of equation, which 
we solve rising a Foiiricr transformation. The then 
obtained vector potential is used (in part, typically 
5%) to update the current vector potential grad- 
ually. The updated vector potential is substituted 
back in the first GL equation and the whole proce- 
dure is repeated until a convergent solution of both 
GL equations is found. 



Temperature dependence 

The temperature dependence of characteristic 
lengths ^ and A and the critical field is incorpo- 



rated in GL theory as ^(T) = 



m 



V|l-T/Teo| 



m 



and Ho2{T) = Hc2{0) 



1- ^ 



X{T) = 
How- 



V|l-T/TeO 

ever, if GL equations are scaled to temperature 
dependent units as given above, that allows us to 
consider any temperature dependence of the pa- 
rameters prior to the start of the simulation. A 
number of works in the past were dedicated to the- 
oretically improved temperature dependence of ^, 
A, and k,. Underlying effects for reported differences 
may be nonlocality, clean/dirty limit effects, strong 
coupling. GL parameter k is temperature indepen- 
dent in GL theory, but for e.g. type-I supercon- 
ductors Ginzburg himself suggested the correction 
k(T) = K(0)/(l-ft2) (t he two -fluid model), Bardeen 
gave k{T) = K(0)/\/l + t^, while Gorkov found 
k{T) = k(0)(1 - 0.24t2 + OMt^), where t = T/T^O. 
Although these formulae go beyond GL theory, we 
emphasize here that they can be implicitly included 
in the calculation, with aim to improve the overall 
validity of the theory. 

It is also worth mentioning that simulations can 
account for the exact experimental procedure with 
respect to cooling. For the zero-field cooled regime 
the initial value of ip should be taken ps 1, and in 
the field-cooled regime opposite, i.e. w since 
superconductivity nucleates from the normal state. 



Fourier transform and 2D approximation 

It is well known that in superconductors the mag- 
netic field penetrates only into a relatively small 
depth A and that screening currents flow in the sur- 
face layer, decaying exponentially in the bulk be- 
yond this length. However, one question arises: what 
happens when the superconductor is thinner than 
the London penetration depth? Tinkham argued, 
as is now generally accepted, that the currents in 
thin superconductors may be considered constant 
over the thickness. Consequently, they have no z- 
component, and the boundary condition (4) is auto- 
matically fulfilled at top and bottom surfaces of the 
sample. Prozorov et al. [14] confirmed Tinkham's 
original assumption, considering in detail the cur- 
rent density (in) homogeneity throughout the thick- 
ness of superconducting films. 
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Therefore, for a superconductor with thickness 
d < A, ^ we assume the uniform distribution of cur- 
rent in the z direction. From the first GL equation 
(2) then follows the same behavior for the order 
parameter, and the 3D problem is reduced to a 2- 
dimensional superconductor. However, an issue of 
solving the GL equation for the vector potential re- 
mains, and here we discuss the details. 

We are actually solving three equations of shape 
— Ac^AA = j. Prom Fourier theory we know that 



oo 

F{k) = J f{x) exp{—2Trixk)dx, 

— oo 
oo 

f{x) = J F{k) exp{2'iTixk)dx. 



(13) 



(14) 



We use this to solve the equation analytically in the 
z-direction. In the x- and y-direction, we solve nu- 
merically using FFT which is based on the following 
Fourier relationship 



Af-l 



JV-1 



fj = ^nsmijn^), 



(15) 



(16) 



where we have chosen the sine transform, to respect 
the boundary condition A{r oo) — >■ 0^ . 

Every component of the vector potential can 
therefore be represented as 

Af-l N-1 

A{x,y,z) = X! m 
i=i j=i 



/ 



dkaij{k) ex.p(2'Kizk) sin(-^) sin(-^), (17) 



where x and y are integer numbers between 1 and 
N and z can be any real number. We assume the 
sample with thickness d to be centered at z = 0. The 
Laplacian of A can be obtained analytically as: 

JV-liV-l °? 

AA{x,y,z)=J2Yl dk{-{2nkf-q^) 

aij{k) exp{2nizk) sin(— ) sin(— ), (18) 



^ The A we solve for here represents solely the field induced 
by the superconductor and does not include the applied field. 



where we introduced q'^ = {iiT/Lx)^+{j7r/Ly)^. L^^y 
is the size of the simulation region in the x- and y- 
direction respectively. 

Now we introduce j{x,y) as the current uniform 

over the sample thickness 

j{x, y, z) = j{x, y)Il{z, -d/2, d/2), 

where 11 represents the step-like function which is 
1 inside the interval [—d/2, d/2] and outside. The 
Fourier transform of 11 yields: 

sin(7rfcrf) 



n(A;) = d- 



Trkd 



Using sine transform, we can express the current as: 

N-l N-1 

3{x,y) = XI H bijsm{iTvx/N)sm{jiTy/N), 

i=l = 1 

with coefficients: 

= iV2 H H 3{x,y)sm{iTT^)sm{jTT^). 

x=\ y=l 

Substituting this into initial GL equation we get: 

,v-i.v-i . . °? 

mi^^'^^'W^^^'^^'^'' I dkexp{2mzk) 

i=l J=l -oo 

K'[{27rkr + q^Mk) - bi.d"^^^^^ = 0, (19) 

which can only be true if the terms in the brackets 
equal zero, i.e. when 

bij sin(7rfcd) 
-J"''^''' = (27rfc)2 + g2 ^kd ■ 

We now revert to the definition of A in terms of its 
Fourier transform, and obtain: 



Af-l 



N-l 



^ A{x,y,z) = ^ sin(i7ra;/iV) ^ sm{j'Ky/N)b^^ 



i=l 



dk exp{2Trizk) 



sm{'jTkd) 



{2'Kkf + g2 -Kkd 



(20) 



The last term can be integrated analytically to get: 
1 



2 (1 — cosh(gz) exp(— dg/2)) , if z < d/2, 



1 

dq 



2 (sinh(dg'/2) exp(-g2;)) , \{ z > d/2. 



We insert 2; = directly in above equation, and 
obtain 5^ (1 - exp(-dg/2)). 
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Recapitulating, we can express A{x, y,z = 0) in 
its Fourier components: 



A{x,y,0) = 

N-l N-1 

^ sin(i7ra;/iV) ^ sm{jTry/N)aij, 
with as coefBcients: 



bij 1 — exp(— fi(7/2) 



(21) 



(22) 



Note that the common definition for efi^ective GL 
parameter for thin superconductors Keff = /d is 
therefore just a lirnMing case obtained for extremely 
thin samples {j{x, y, z) — S{z)i{x, y)). 

Alternatively, one can calculate the mean vector 
potential (averaged in ^-direction) inside the sam- 
ple, instead of the values in z = plane only. In that 
case, the Fourier coefficients become: 



1 bin /, „sinh(dg/2)\ , . ,„x 
^ ' d- 2 exp{-dij/2). 



The anisotropy 

Modern fabrication methods enable experiments 
on various hybrid structures. Among others, one can 
think of complex 3D superconductor-metal hybrids, 
as well as samples made of different superconduct- 
ing materials. In addition, materials such as high- 
Tc superconductors are mostly layered, and exhibit 
clear anisotropy in different directions. In what fol- 
lows, we show how those are incorporated in our GL 
formalism. 

To begin with, we consider a superconductor- 
metal hybrid. As already pointed out by deGennes, 
the leakage of Cooper-pairs from the superconduc- 
tor into the metal can be modeled through the 
modified boundary condition for a superconductor- 
normal metal interface as 



n ■ {-ifN - -^^)'>P 



(23) 



interface 



where quantity b is positive and measures the dis- 
tance outside the boundary (in the normal metal) 
where the order parameter becomes zero if the slope 
at the interface is maintained. As a consequence, 
the discrcitizcd first GL equation (10) must be mod- 
ified at the boundary of the superconductor, and 
the terms 'tpj/af^ perpendicular to the interface get 
a prefactor (1 — 6). 



In the case of a superconductor-superconductor 

hybrid, the situation is far more complex. The sim- 
plest scenario is that only is different between two 
materials, in which case the only parameter a varies 
in the sample. This is put in GL equations directly, 
and generates no extra terms in the discretization. 
Note that the crossing currents are also calculated 
directly, and no special boundary condition (such as 
Eq. (23)) is needed. In the case of varied mean free 
path, both a and (3 change (proportional to I and 
P respectively), but this still docs not generate ad- 
ditional terms in GL equations and can be put in 
the calculations directly. However, to consider mass 
anisotropy in the sample, one must go back to the 
very derivation of GL equations to establish that 

an additional term V— — ^ ( — zV — A] ip appears 

'^\'^ -iV 1^ ) y J 

on the right side of Eq. (2). This term is then dis- 
cretized using Eq. (12) and added to Eq. (10) in the 
numerical procedure. 

Note that for real layered samples (such as high- 
Tc ones) the coupling between the layers must also 
be introduced into GL equations. This is the known 
Lawrence-Doniach extension to the GL model, 
where Josephson coupling is active between the su- 
perconducting layers separated by an insulator of 
thickness rfj. As a result, at planes z = ±di/2 a term 
VP±dj2 - V'T<i,/2exp(±iAf„J]/d, is added to the 

first GL equation (here A^^^ = J^^^^^^^A^dz). The 
second GL equation remains the same, but an ad- 
ditional equation describing the Josephson current 
enters the calculation at all S-I boundaries as 



2di 



di/2 



.(24) 



Two-band superconductivity 



Multi-band superconductivity has stirred great 
interest in the scientific community following the 
discovery and further studies of MgB2 [15], having 
the highest critical temperature for a non-copper- 
oxide and non-fullerene superconductor. This mate- 
rial and some others such as boro-carbides can be 
described by two superconducting order parameters 
(hence the name 'two-band' or 'two-gap' supercon- 
ductor). However the real boom in two- and multi- 
band superconductivity followed the last year's dis- 
covery of iron pnictides [16]. The concentration of 
different dopants and applying of external pressure 
can tune the electronic, magnetic and structural 
properties of these materials. Provided that their 
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critical temperature is increased above liquid air, the 
impact of these materials on technology is already 
seen as a 'new iron age'. 

In cases when microscopic structure of the ma- 
terial is not relevant, the Ginzburg-Landau theory 
can describe two-gap superconductors. The first GL 
equation, for each of the two order parameters now 
reads: 

i-iV - Af^i - (xi - |Vi|')V'i - |V'2 

(-iV- a)%2 = (25) 

1 s 

-(-iV - Aftp2 - (X2 - |V'2n tp2 - — Vl 

a ma 



(-zV -Aj ^1 = 0, (26) 



where Xi = 1 — T /T^ {T^ being the critical temper- 
ature of gap i). a = {£,w/i2of and 5 = i/'io/V'2o = 



/9ia20 



, all measured at T = and in absence of 

magnetic field and coupling. 7 and 77 are the cou- 
pling constants representing Josephson and drag ef- 
fect respectively, m = mi /m2 is the ratio of the ef- 
fective masses in two gaps. Lengths are expressed in 
^10, and the order parameters in tpiQ. 

The second GL equation is a single equation, with 

contribution from both condensates. We write it in 

2 

the following form using the relation ^ = ^ : 



-AA = 



1 



^7^ 



^ 1 

■rixi n 



+ T1 



i,l - Aj Vi 

\p* (_iv - a) V2' 

ri (-iV - i) ^2 

-n V2 (-*v - A) V'l 



m K1K2 

Tx 1 



m K1K2 



(27) 



where we singled out the GL parameters for both 
condensates, although this is not realistic. However, 
the latter equation makes directly visible that the 
magnetic behavior of the sample will depend not 
only on ni and K2, but also on the ratio of coherence 
lengths in two condensates (parameter a) , and this is 
of direct relevance to the recently observed type- 1.5 
superconductivity [17]. Note also that all terms in 
Eqs. (25-27) were already present in the same shape 
in the single-gap GL equations. Therefore, Eqs. (25- 
27) can be discretized and solved using the same 
numerical scheme as given in preceding sections. 



In summary, we have explained in more detail 

the numerical approach that discrctizcs Ginzburg- 
Landau equations on a three dimensional Cartesian 
grid. We further reviewed its applications on su- 
perconducting samples of different topologies, the 
hybrid systems, anisotropic samples, layered super- 
conductors, and two-gap superconductors. For thin 
samples, we show the modification of the approach 
that allows fast solving for order parameter in the 2D 
plane, while the vector potential is generally solved 
in all three dimensions. We therefore hope to have 
presented a brief manual which will be of use for 
young researchers in the years to come. 
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